This report runs base, AR1, and AR2 models (variable ~ time) for annual SEV meteorological data using the nlme R package. Mean annual air temperature, minimum annual air temperature, maximum annual air temperature, and total annual precipitation are modeled. More sophisticated models will likely be required after this preliminary analysis.


library(tidyverse)
library(lubridate)
library(nlme)
library(MuMIn)
library(gridExtra)
library(kableExtra)

# Load yearly data ---------------------------------------------

path_to_files <- "../data/processed_data/"

# load and only keep vars of interest
y <- read_csv(paste0(path_to_files, "met_yearly_gap_filled.csv")) %>% 
  mutate(sta = as.factor(sta)) %>% 
  select(sta:ppt)
str(y)
## tibble [110 × 6] (S3: tbl_df/tbl/data.frame)
##  $ sta   : Factor w/ 4 levels "40","42","49",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ year  : num [1:110] 1988 1989 1990 1991 1992 ...
##  $ airt  : num [1:110] 13.8 13.8 13.8 12.5 12.3 ...
##  $ maxair: num [1:110] 35 35 35 37.7 38.6 ...
##  $ minair: num [1:110] -7.94 -7.94 -7.94 -13.76 -15.11 ...
##  $ ppt   : num [1:110] 151 140 232 335 241 ...
summary(y)
##  sta          year           airt           maxair          minair       
##  40:34   Min.   :1988   Min.   :11.91   Min.   :32.01   Min.   :-29.930  
##  42:33   1st Qu.:2001   1st Qu.:13.15   1st Qu.:37.13   1st Qu.:-16.495  
##  49:23   Median :2008   Median :13.82   Median :38.70   Median :-14.780  
##  50:20   Mean   :2007   Mean   :13.89   Mean   :38.50   Mean   :-14.900  
##          3rd Qu.:2015   3rd Qu.:14.67   3rd Qu.:39.84   3rd Qu.:-12.982  
##          Max.   :2021   Max.   :16.43   Max.   :48.72   Max.   : -3.797  
##       ppt       
##  Min.   :122.0  
##  1st Qu.:220.9  
##  Median :247.9  
##  Mean   :262.5  
##  3rd Qu.:301.7  
##  Max.   :599.4

Prepare data for modeling -

Add a time variable for nlme modeling and split data into individual met stations.

# function to select station and variable and add 'time' var
prepare_data_for_nlme <- function(data, sta) {
  # data = dataframe to use
  # sta = met station (in quotes)
  data %>% 
    filter(sta == {{ sta }}) %>% 
    mutate(time = as.numeric(as.factor(year)))
}


m40 <- prepare_data_for_nlme(y, "40")
m42 <- prepare_data_for_nlme(y, "42")
m49 <- prepare_data_for_nlme(y, "49")
m50 <- prepare_data_for_nlme(y, "50")

Set up functions -

This code produces several functions that will be useful for producing cleaner code.

# function for plotting preliminary yearly graph -----------------------------------------------------------

plot_yearly_prelim <- function(data, var, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = year, y = {{ var }}, col = sta)) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    facet_wrap(~ sta) +
    labs(title = title,
         x = "Year",
         y = y_axis_label)
}

# function for plotting graphs with best nlme model results ------------------------------------------------

plot_yearly_results <- function(data, var, coeffs, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # coeffs = object containing coefficients and p-values of best model
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = year, y = {{ var }})) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    labs(title = title,
         x = "Year",
         y = y_axis_label,
         caption = paste("slope = ", round(coeffs[2,1], 3), "\n(p-value = ", coeffs[2,2], ")"))
}


# functions to run nlme models ----------------------------------------------------------

# base model function
base_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      method = "ML",
      na.action = na.omit)
}



# AR1 model function
AR1_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 1),
      method = "ML",
      na.action = na.omit)
}


# AR2 model function
AR2_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 2),
      method = "ML",
      na.action = na.omit)
}

Mean Annual Air Temperature (C) -

plot_yearly_prelim(y, airt, "Annual Mean Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "airt")
m40_mAR1 <- AR1_nlme_model(m40, "airt")
m40_mAR2 <- AR2_nlme_model(m40, "airt")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_mbase 3 53.68867
m40_mAR1 4 55.62003
m40_mAR2 5 58.37968
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   52.88867 57.46775 -23.44433
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 13.271801 0.17431449 76.13711  0.0000
## time         0.021443 0.00868861  2.46794  0.0191
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1560729 -0.5446061  0.2023188  0.6881964  1.6138251 
## 
## Residual standard error: 0.4821986 
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, airt, m40_coeffs, "Met 40 - Deep Well \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "airt")
m42_mAR1 <- AR1_nlme_model(m42, "airt")
m42_mAR2 <- AR2_nlme_model(m42, "airt")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_mbase 3 50.38362
m42_mAR1 4 52.94158
m42_mAR2 5 55.65049
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   49.55604 54.04556 -21.77802
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 12.618607 0.17205515 73.34048  0.0000
## time         0.011535 0.00883012  1.30629  0.2011
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -1.987976067 -0.644320891 -0.002985058  0.916472339  1.576226874 
## 
## Residual standard error: 0.468135 
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, airt, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "airt")
m49_mAR1 <- AR1_nlme_model(m49, "airt")
m49_mAR2 <- AR2_nlme_model(m49, "airt")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_mbase 3 34.60607
m49_mAR1 4 37.16381
m49_mAR2 5 39.81788
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   33.34291 36.74939 -13.67145
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 14.203751 0.19776763 71.82041  0.0000
## time         0.024676 0.01442369  1.71081  0.1018
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1832629 -0.6706411  0.2109746  0.6989873  1.7965345 
## 
## Residual standard error: 0.4384421 
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, airt, m49_coeffs, "Met 49 - Five Points \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "airt")
m50_mAR1 <- AR1_nlme_model(m50, "airt")
m50_mAR2 <- AR2_nlme_model(m50, "airt")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_mbase 3 33.37573
m50_mAR1 4 35.87783
m50_mAR2 5 39.37225
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   31.87573 34.86293 -12.93786
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 14.831383 0.22625535 65.55152  0.0000
## time         0.051514 0.01888743  2.72741  0.0138
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81172625 -0.64841217 -0.09216346  0.63402578  1.96854695 
## 
## Residual standard error: 0.462067 
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, airt, m50_coeffs, "Met 50 - Blue Grama \nAnnual Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Mean Annual Air Temperature

# airt results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'


Minimum Annual Air Temperature (C) -

plot_yearly_prelim(y, minair, "Annual Minimum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "minair")
m40_mAR1 <- AR1_nlme_model(m40, "minair")
m40_mAR2 <- AR2_nlme_model(m40, "minair")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_mbase 3 189.6589
m40_mAR1 4 191.1119
m40_mAR2 5 193.7564
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   188.8589 193.438 -91.42944
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept) -13.656159 1.2874555 -10.607092  0.0000
## time         -0.136817 0.0641725  -2.132025  0.0408
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6474716 -0.3773667  0.0234212  0.4512047  1.7215955 
## 
## Residual standard error: 3.561432 
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, minair, m40_coeffs, "Met 40 - Deep Well \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "minair")
m42_mAR1 <- AR1_nlme_model(m42, "minair")
m42_mAR2 <- AR2_nlme_model(m42, "minair")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_mbase 3 189.4901
m42_mAR1 4 190.6799
m42_mAR2 5 192.9317
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   188.6625 193.152 -91.33124
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept) -13.220095 1.4158546 -9.337184  0.0000
## time         -0.094695 0.0726637 -1.303196  0.2021
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.52824489 -0.32599081 -0.05008804  0.56078847  2.47065756 
## 
## Residual standard error: 3.852318 
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, minair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "minair")
m49_mAR1 <- AR1_nlme_model(m49, "minair")
m49_mAR2 <- AR2_nlme_model(m49, "minair")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_mbase 3 123.2896
m49_mAR1 4 125.6878
m49_mAR2 5 128.9351
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   122.0264 125.4329 -58.01321
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept) -14.590316 1.3596668 -10.730802    0.00
## time          0.006294 0.0991639   0.063475    0.95
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.0179895 -0.3267865  0.1783908  0.5709077  1.3371903 
## 
## Residual standard error: 3.014321 
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, minair, m49_coeffs, "Met 49 - Five Points \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "minair")
m50_mAR1 <- AR1_nlme_model(m50, "minair")
m50_mAR2 <- AR2_nlme_model(m50, "minair")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_mbase 3 112.1543
m50_mAR1 4 115.2297
m50_mAR2 5 118.1863
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   110.6543 113.6415 -52.32717
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept) -14.264158 1.6215367 -8.796691  0.0000
## time          0.072444 0.1353633  0.535179  0.5991
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8049348 -0.3200412  0.1852600  0.5306493  1.1520424 
## 
## Residual standard error: 3.311562 
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, minair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Minimum Annual Air Temperature

# minair results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'


Maximum Annual Air Temperature (C) -

plot_yearly_prelim(y, maxair, "Annual Maximum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "maxair")
m40_mAR1 <- AR1_nlme_model(m40, "maxair")
m40_mAR2 <- AR2_nlme_model(m40, "maxair")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_mbase 3 140.3235
m40_mAR1 4 135.6487
m40_mAR2 5 136.8874

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m40_mAR1)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   134.2694 140.3748 -63.13468
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##       Phi 
## 0.4530549 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 37.64031 0.9617813 39.13604  0.0000
## time         0.08334 0.0473790  1.75903  0.0881
## 
##  Correlation: 
##      (Intr)
## time -0.862
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6634576 -0.5746671 -0.1493423  0.6655313  3.1746966 
## 
## Residual standard error: 1.732308 
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mAR1)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, maxair, m40_coeffs, "Met 40 - Deep Well \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "maxair")
m42_mAR1 <- AR1_nlme_model(m42, "maxair")
m42_mAR2 <- AR2_nlme_model(m42, "maxair")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_mbase 3 129.8644
m42_mAR1 4 127.7573
m42_mAR2 5 130.0821

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m42_mAR1)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   126.3287 132.3147 -59.16435
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##       Phi 
## 0.3867261 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 34.94019 0.8290183 42.14647  0.0000
## time         0.06777 0.0421570  1.60763  0.1181
## 
##  Correlation: 
##      (Intr)
## time -0.864
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.94732099 -0.54279907 -0.08489296  0.66776673  2.67346131 
## 
## Residual standard error: 1.572183 
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mAR1)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, maxair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "maxair")
m49_mAR1 <- AR1_nlme_model(m49, "maxair")
m49_mAR2 <- AR2_nlme_model(m49, "maxair")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_mbase 3 116.1389
m49_mAR1 4 117.1260
m49_mAR2 5 116.5718

The base model has the lowest AICc and is used for plotting.

summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   114.8758 118.2822 -54.43788
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 41.07945 1.1639166 35.29415  0.0000
## time        -0.10151 0.0848873 -1.19584  0.2451
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3311794 -0.5068364 -0.1698531  0.1037714  3.2757736 
## 
## Residual standard error: 2.580352 
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, maxair, m49_coeffs, "Met 49 - Five Points \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "maxair")
m50_mAR1 <- AR1_nlme_model(m50, "maxair")
m50_mAR2 <- AR2_nlme_model(m50, "maxair")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_mbase 3 63.96627
m50_mAR1 4 63.49305
m50_mAR2 5 66.10862

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m50_mAR1)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   60.82638 64.80931 -26.41319
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##        Phi 
## -0.4018488 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 38.90612 0.3271151 118.93708  0.0000
## time         0.07382 0.0274879   2.68568  0.0151
## 
##  Correlation: 
##      (Intr)
## time -0.882
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.65156196 -0.77367285 -0.05306683  0.49518631  2.02742794 
## 
## Residual standard error: 0.9854843 
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mAR1)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, maxair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Maximum Annual Air Temperature

# maxair results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'


Total Annual Precipitation (mm) -

plot_yearly_prelim(y, ppt, "Total Annual Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "ppt")
m40_mAR1 <- AR1_nlme_model(m40, "ppt")
m40_mAR2 <- AR2_nlme_model(m40, "ppt")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m40_mbase 3 376.9800
m40_mAR1 4 379.0143
m40_mAR2 5 380.8387
summary(m40_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##      AIC      BIC  logLik
##   376.18 380.7591 -185.09
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 239.73384 20.234522 11.847764  0.0000
## time         -0.79358  1.008579 -0.786827  0.4372
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.7633224 -0.6843958  0.0941523  0.8117423  1.7513612 
## 
## Residual standard error: 55.97387 
## Degrees of freedom: 34 total; 32 residual
m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]
(m40_plot <- plot_yearly_results(m40, ppt, m40_coeffs, "Met 40 - Deep Well \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "ppt")
m42_mAR1 <- AR1_nlme_model(m42, "ppt")
m42_mAR2 <- AR2_nlme_model(m42, "ppt")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m42_mbase 3 386.6521
m42_mAR1 4 389.0729
m42_mAR2 5 390.7934
summary(m42_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   385.8245 390.314 -189.9122
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 391.6703 28.079405 13.948670  0.0000
## time         -3.5405  1.441075 -2.456826  0.0198
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.60096602 -0.53817435 -0.06228751  0.56988895  3.13557611 
## 
## Residual standard error: 76.39964 
## Degrees of freedom: 33 total; 31 residual
m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]
(m42_plot <- plot_yearly_results(m42, ppt, m42_coeffs, "Met 42 - Cerro Montoso \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "ppt")
m49_mAR1 <- AR1_nlme_model(m49, "ppt")
m49_mAR2 <- AR2_nlme_model(m49, "ppt")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m49_mbase 3 253.8175
m49_mAR1 4 256.5620
m49_mAR2 5 259.8661
summary(m49_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   252.5543 255.9608 -123.2772
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 254.18442 23.215104 10.949096  0.0000
## time         -1.99381  1.693136 -1.177584  0.2521
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.01258367 -0.67241549  0.07653863  0.62248810  2.09072501 
## 
## Residual standard error: 51.46686 
## Degrees of freedom: 23 total; 21 residual
m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]
(m49_plot <- plot_yearly_results(m49, ppt, m49_coeffs, "Met 49 - Five Points \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "ppt")
m50_mAR1 <- AR1_nlme_model(m50, "ppt")
m50_mAR2 <- AR2_nlme_model(m50, "ppt")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")
df AICc
m50_mbase 3 222.6721
m50_mAR1 4 225.6402
m50_mAR2 5 229.2037
summary(m50_mbase)
## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC   logLik
##   221.1721 224.1593 -107.586
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 245.26421 25.695541 9.545011  0.0000
## time          0.26436  2.145022 0.123244  0.9033
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.95082391 -0.49761448 -0.02635629  0.88260113  1.46242374 
## 
## Residual standard error: 52.47638 
## Degrees of freedom: 20 total; 18 residual
m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]
(m50_plot <- plot_yearly_results(m50, ppt, m50_coeffs, "Met 50 - Blue Grama \nTotal Annual Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Total Annual Precipitation

# ppt results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'